#include "error.cuh"
#include <stdio.h>

#ifdef USE_DP
    typedef double real;
#else
    typedef float real;
#endif

const int NUM_REPEATS = 100;
const int N = 100000000;
const int M = sizeof(real) * N;
const int BLOCK_SIZE = 128;


void timing(real *h_x, real *d_x, const int method);


void __global__ reduce_global(real *d_x, real *d_y){
    int tid=threadIdx.x;
    real* x=blockIdx.x*BLOCK_SIZE+d_x;
    for(int i=BLOCK_SIZE>>1;i>0;i>>=1){
        if(tid<i)
            x[tid]+=x[tid+i];
        __syncthreads();
    }

    if(tid==0)
        d_y[blockIdx.x]=x[0];

}





void __global__ reduce_shared(real *d_x, real *d_y){
    const int tid=threadIdx.x;
    const int idx=threadIdx.x+blockIdx.x*BLOCK_SIZE;
    __shared__ real s_y[128];
    s_y[tid]=(idx<N)?d_x[idx]:0.0;
    __syncthreads();
    for(int offset=BLOCK_SIZE>>1;offset>0;offset>>=1){
        if(tid<offset){
            s_y[tid]+=s_y[tid+offset];
        }
        __syncthreads();
    }
    if(tid==0)
        d_y[blockIdx.x]=s_y[0];
        //atomicAdd(d_y,s_y[0]);
};


void __global__ reduce_dynamic(real *d_x, real *d_y){

    const int tid=threadIdx.x;
    const int idx=threadIdx.x+blockIdx.x*BLOCK_SIZE;
   
    extern __shared__ real s_y[];
    s_y[tid]=(idx<N)?d_x[idx]:0.0;
    __syncthreads();
    for(int offset=BLOCK_SIZE>>1;offset>0;offset>>=1){
        if(tid<offset){
            s_y[tid]+=s_y[tid+offset];
        }
        __syncthreads();
    }
    if(tid==0)
        d_y[blockIdx.x]=s_y[0];
};

real reduce(real *d_x, const int method){
    int grid_size=(N+BLOCK_SIZE-1)/BLOCK_SIZE;
    const int ymem=grid_size*sizeof(real);
    const int smem=BLOCK_SIZE*sizeof(real);
    real* d_y;
    CHECK(cudaMalloc(&d_y, ymem));
    real* h_y=(real*)malloc(ymem);
    switch (method)
{
    case 0:
        reduce_global<<<grid_size, BLOCK_SIZE>>>(d_x, d_y);
        break;
    case 1:
        reduce_shared<<<grid_size, BLOCK_SIZE>>>(d_x, d_y);
        break;
    case 2:
        reduce_dynamic<<<grid_size, BLOCK_SIZE, smem>>>(d_x, d_y);
        break;
    default:
        printf("Error: wrong method\n");
        exit(1);
        break;
}
    CHECK(cudaMemcpy(h_y, d_y, ymem, cudaMemcpyDeviceToHost));
    real result=0.0;
    for(int i=0; i<grid_size; i++)
        result+=h_y[i];
    free(h_y);
    CHECK(cudaFree(d_y));
    return result;
}


int main(void)
{
    real *h_x = (real *) malloc(M);
    for (int n = 0; n < N; ++n)
    {
        h_x[n] = 1.23;
    }
    real *d_x;
    CHECK(cudaMalloc(&d_x, M));

    printf("\nUsing global memory only:\n");
    timing(h_x, d_x, 0);
    printf("\nUsing static shared memory:\n");
    timing(h_x, d_x, 1);
    printf("\nUsing dynamic shared memory:\n");
    timing(h_x, d_x, 2);

    free(h_x);
    CHECK(cudaFree(d_x));
    return 0;
}

void timing(real *h_x, real *d_x, const int method){
    real sum=0;
    for(int i=0;i<NUM_REPEATS;i++){
        CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));

        cudaEvent_t start, stop;
        CHECK(cudaEventCreate(&start));
        CHECK(cudaEventCreate(&stop));
        CHECK(cudaEventRecord(start));
        cudaEventQuery(start);

        sum = reduce(d_x, method);

        CHECK(cudaEventRecord(stop));
        CHECK(cudaEventSynchronize(stop));
        float elapsed_time;
        CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));
        printf("Time = %g ms.\n", elapsed_time);

        CHECK(cudaEventDestroy(start));
        CHECK(cudaEventDestroy(stop));

    }
    printf("sum = %f.\n", sum);
}

